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£Sj ■ We present an exact mathematical description of the Meissner effect and of the vortex state in 

' periodic thin-layer superconductor/insulator structures with an arbitrary number of identical junc- 

tions N — 1 (2 < TV < oo, where N is the number of superconducting layers) in the presence of a 
rjr_^ ■ static parallel external field H. Based on an analytical analysis of the coupled static sine-Gordon 

(SG) equations for the phase differences, we obtain a complete classification of all possible types 
of physical solutions. We prove that at H > these equations admit only Meissner solutions and 
topological "vortex-plane" solutions. Both the types of solutions are characterized, in general, by 
[t] Josephson lengths ([-y] is the integer part of -y). We derive an explicit analytical expres- 
sion for the superheating field of the Meissner state, H s , as a function of N and show that H„ 
simultaneously determines the penetration field for the vortex planes. For H <C H s , we obtain a 
^ , closed-form analytical expression for the Meissner field and investigate a transition to the infinite 

S_h ■ (N = oo) layered-superconductor limit. Thermodynamically stable "vortex-plane" solutions repre- 

Oh' sent coherent chains of Josephson vortices (one vortex per each insulating layer in a chain). Being a 

natural generalization of ordinary Josephson vortices in a single junction, the vortex-plane solutions 
inherit such properties of the former as periodicity along the layers and the overlapping of states 
with different topological numbers. We obtain exact analytical expressions for the self-energy of a 
vortex plane and for the lower critical field H c \. In contrast to a prevailing view, the coupled SG 
equations do not possess any single-vortex solutions for H > 0, as well as such vortex solutions as 
[ the "triangular lattice" and the "row mode". Thus, single- vortex solutions appear only in the limit 

H = 0, for which case we provide a detailed analytical description. The general consideration is 
illustrated by two exactly-solvable examples (iV = 2,3). Experimental implications of the results 
are discussed. 

PACS numbers: 74.50. +r, 74.80.Dm 
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We present a rigorous mathematical examination of the problem of the Meissner effect and of the vortex structure 
in thin-layer Josephson-junction stacks and layered superconductors, with a static, homogeneous external magnetic 
field H > applied parallel to the layers (along the z axis, see Fig. 1.) We consider periodic systems composed of an 
arbitrary number N — 1 of identical superconductor/insulator (S/I) junctions (2 < N < oo, where N is the number 

S of S-layers, with x being the layering axis), occupying an arbitrary region [— L, L] along the y axis. 
We begin with the investigation of analytical properties of a finite set of coupled static sine-Gordon (SG) equations 
for phase differences </>„(l < n < N—l), obtained by the minimization of a microscopic Gibbs free-energy functionaloBI 
Equations of this type were first derived within the framework of a phenomenological approach.B Undet-|Corresponding 
redefinition of the parameters, they also apply to the phenomenological Lawrence-Doniach (LD) made] 3 with N < oo. 
Although particular numerical solutions to these equations were obtained in several publications,™ - !^ any analytical 
. analysis of their general properties has not been undertaken up until now. r-. 

Making use of standard methods of the theory of differential equations,El we arrive at a complete classification of 
the solutions to the SG equations subject to appropriate physical boundary conditions at y = ±L. In particula: 
we prove that for H > these equations admit only Meissner solutions and topological "vortex-plane" solutions a 
Both the types of solutions are characterized, in general, by [y] Josephson lengths \jk ([y] ^ s the integer part of ^ 
k = 0, 1, . . . , [^1 — I). The Meissner solutions persist up to a certain superheating field of the Meissner state, H s l- We 
show that the field H s l simultaneously determines the penetration field for the vortex planes. For L 3> Aj max = Ajo, 
we derive an explicit analytical expression for the superheating (penetration) field H s = H soo as a function of N > 2. 
For H <§; H s , we obtain a closed-form analytical expression for the Meissner local field H n (y), valid for any N > 2, 
and investigate a transition to the infinite (N = oo) layered-superconductor limit. 

Thermodynamically stable vortex-plane solutions, physically, represent coherent chains of Josephson vortices (one 
vortex per each I-layer in a chain), positioned in planes parallel to the coordinate axes x, z. (See Fig. 2.) Such solutions 
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were thoroughly studied for infinite layered superconductors in Refs. . For H c \ <§C H <§; H C 2, thew correspond to 
the "transparent state" considered by Theodorakis within the framework of the infinj^UXrnodel.EI In the case of 
a double-junction stack (N = 3), they coincide with the well-known "in-phase mode" .uli3 tlHa Significantly, coherent 
Josephson vortex configurations that can be identified with the vortex-plane solutions have been directly observed 
in experiments on artificial double-junction stacksE^I and weakly-coupled multilayers Besides giving a proof of the 
existence and stability of the vortex-plane solutions for 3 <iV < oo, we discuss their main physical properties, such as 
periodicity in the y direction and the overlapping of states with different numbers of vortex planes N v . In the limit 
L 3> A,/ max = A jo, we derive exact analytical expressions for the self-energy E v of the state with N v = 1 and for the 
lower critical field H c \. We also investigate a physical and mathematical relationship between vortex-plane solutions 
for N > 3 and the ordinary Josephson vortices in a single junction (N — 2). 

In contrast to a prevailing view, the coupled SG equations do not passes any single- vortex solutions for H > and 
L < oo, as well as such vortex solutions as the "triangular lattice"EJ~t3 and the "row mode" .19 Thus, single- vortex 
solutions appear only in the limit H = 0, L = oo, which, physically, signifies their absolute thermodynamic instability. 
In this regard, it should be emphasized that single-vortex configurations were originally pxnpnsed-.for infinito-Lawped 
superconductors without an appropriate analysis of the actual conditions of their existence. LDii3ll3'E3 (See Refs.E2rEJl3 for 
the criticism.) An excessive preoccupation with single- vortex configurations in previous publications can be explained 
by the confusion of the problem of finding the lowest-energy topological solutions at H — with the problem of the 
minimization of the Gibbs free energy at H > 0: Although at H = 0, L = oo the energy of single- vortex configurations 
is lower than that of the vortex planes, the former ones are irrelevant to the problem of minimization because of their 
absolute instability at H > 0, L < oo. For the same reason, single-vortex solutions for H = 0, L = oo cannot be used 
for estimates of the lower critical field H cl . The absence of single- vortex configurations at H > was established 
for infinite layered superconductors by the use of exact variational methods in Refs.BH . In order to clarify this issue 
for N < oo, we analyze those features of single- vortex solutions at H = 0, L = oo which preclude their existence at 
H > 0, L < oo. 

Section II of the paper is devoted to exact mathematical formulation of the problem. In section III, we derive 
all major mathematical and physical results sketched above. Mathematical Lemmas 1, 2 and the Theorem of 
subsection IIIA form the basis for a subsequent physical analysis in subsections IIIB-IIID. The general consideration 
of section III is illustrated in section IV by two exactly-solvable examples, namely of a single thin-layer junction 
(JV = 2) and of a double-junction stack (N = 3). (In section IV, we present a new exact solution to the single SG 
equation, valid for arbitrary H > and L > 0.) The obtained results are discussed in section V. Appendices A-C 
contain mathematical derivations and proofs omitted in the main text. 



II. FORMULATION OF THE PROBLEM 



We consider a periodic structure consisting of alternating N superconducting(S) and N — 1 insulating (I) layers 
(2 < N < oo) in a parallel, static, homogeneous external magnetic field < H <C H C 2, where H C 2 is the upper critical 
field. (See Fig. 1.) The S-layer thickness is a, and p is the period; the length of the structure in the y direction is 
W = 2L, and W z is the length of the structure in the z direction (W z — * oo). We set % = c = 1 and assume that 

Tc0 ~ T «l, (2.1) 



Co « a, (2.2) 



a <C mm 



{C(T),A(T),a-%}, 



3tt 2 
7C(3) 



dttD(t) < 1, 



(2.3) 



a«j). 



(2.4) 



Here T c0 is the critical temperature of an isolated S-layer, £ is the BCS coherence length; ((T) and X(T) are 
respectively, the Ginzburg-Landau (GL) coherence length and the penetration depth; D(cos8) is the incidence- angle- 
dependent tun neli ng pro bab ility of the I-layer between two successive S-layers 



Conditions (2.1) and (2.2) ensure the applicability of the GL-type expansion within each S-layer. Condition (2.3) 
corresponds to the thin S-layer limit, whereas condition (2.4) is employed here for the sake of mathematical simplicity 
only. To simplify further the mathematical description, wc introduce dimensionless units by 
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where the quantities on the left-hand side are dimensional, with Aj = (87rejop) being the Josephson penetration 
depth (Jo is the density of the Josephson current in a single junction with thick electrodes) and H s — (epAj) being 
the superheating (penetration) field of the infinite layered superconductor oB In our dimcnsioiilees units, for example, 
the flux quantum is $ = ft, and the lower critical field of the infinite layered superconductorEJB is H c \ = — . 

Under the above conditions, the structure is completely described by a closed set of self-consistent microscopic 
equations for the reduced modulus of the superconducting order parameter in the nth S-layer f n (y) (0 < f n < 1, 
n = 0, 1, . . . ,N — lX-and the phase difference between the nth and (n — l)th S-layers 4> n {y) {<frn = ~ fn-ii 
n = 1,2,..., N- 1)M 



fo(y) - / 3 (y) - r(T) 



e 2 d 2 f (y) , 2[H-H l {y)Y , 1 



2 dy 2 



t 2 f$(y) 



2 [fo(y) ~ fi (v) cos (pi (y)} 



fn(y) - f n {y) = r(T) 



e 2 d 2 f n (y) 2[H n (y)-H n+1 (y)Y 
2 dy 2 + e*fl{y) 



"2 Pfn(y) - fn+i(y) cos0„ +1 (y) - f n -i(y) cosc/> n (y)} 



, 1 < n < N - 2, 



(2.5) 



f N -i(y) - fl-^y) = r(T) 



e 2 d 2 f N _ 1 (y) 2[H-H N ^(y)Y 
2 dy 2 + e 2 /^!(y) 



+ 7T [/jV-l(y) ~ /jV-2(y) COS0Ar_i(y)] 



-^(±£)=0, 0<n<7V-l; 



(2.6) 



- H„(j/)] - -j^-j^ [H n (y) - H n _i(y)] - e 2 ff„(y) 



/„ 2 _i(y) 



e 2 #n(^) 
2 rfy : 



1 < n < TV - 1, 



(2.7) 
(2.8) 



dy 



- (±L) = 2H, 1 < n < N — 1, 



(2.9) 



where 
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r(T) 



C 2 (T)a 



rip 



and the local magnetic field in the nth I-layer (n — 1 < x < n) is given by 

1 V 

H n {y) = 2 du fn{u)fn-i(u) sin (f) n (u) + H 



duf n (u)f n -i(u) sin<p n (u) + H. 



(2.10) 



Note that although in most physical situations 6 « 1, we will not nee d the smallness of e in our mathematical 
consideration. Because of the obvious property f n (y) = fn(~V)> equation ( 2.1C ) implies that the phase differences <j> n 
meet the condition 



4>n(y) = -4>n(-y) + 0mod27r. 

The dimensionless Gibbs free energy £1(11), normalized via the relation 

AttQ,(H) 
H2(T)a\ Joo W z mh 

in terms of the mean- held quantities f n , <p n and H n (y) has the form 

N-l L 



(2.11) 



- [H n+1 (y) H n {y)f + 2r(T) [H n (y) H} 2 



r(T) 



N-l H. 

E J dy[f 2 n _ 1 {y) + fl{y)-V n{y)f n -i{y) cos <j> n {y)\ 



(2.12) 



Here, the sum of the three terms in the first line on the right-hand side is the condensation energy, the sum of the 
two terms in the second line is the electromagnetic energy E em , and the last term is the Josephson energy Ej. The 
intralayer current in the nth S-layer J n (y) (normalized to H s ) and the density of the Josephson current between the 
nth and the (n — l)th S-layers j n ,n-i(y) (normalized to jo) are given by 



J n (y) = -L [H n (y) - H n+1 (y)} , < n < N - 1, 

47T 



and 



3n,n-l(y) = 2^fM = / n (y)/ n _ 1 (y) sin0„(y), 1 < n < N - 1, 
dy 

respectively. 

Assuming that the temperature range satisfies the condition of the weak-coupling limit 

r(T) « 1, 



(2.13) 



(2.14) 



(2.15) 



one can obt ain a perturbative solution for /„ and <p n up to any desired order in r(T), starting from the zero-order 
solution to (|J), (|J), 



4 



fn = 1, 



(2.16) 



and the zero-order equations for <fi n , 

H n+1 (y) - (2 + e 2 ) H n (y) + H n ^(y) 



e 2 d(f> n {y) 
2 dy l 



1< n< N- 1, 



(2.17) 



where H n (y) are given by (2.10) with /„ = 1 and satisfy the boundary conditions (2.8). For most applications, it is 
s uffic ient to consider express ions for phy sical quantities only in leading order in r(T). Thus, for example, substituting 
( [2.16 ) and the solution of ( [2.17 ) into (2.12) immediately yields a first-order expansion for the Gibbs free energy, 
because first-order corrections to the condensation-energy term cancel out. 

A detailed mathematical analysis of Eqs. ( [2.17 ) is the subject of section III. Here we point out that these equations 
can be transformed into a very useful for application form by solving for H n (y) (see Appendix A for mathematical 
details): 



H n (y) = h n {y) + H n , 



(2.18) 



K{y) = G ( n > m ) dy 



(2.19) 



,-N+n 



N-T 



,-N 



(2.20) 



where G(n,m) are given by (A9), and fi is given by (|A5|). By (|A8|), and (2J)), ( |A13| ), expression ( |2.1g) exp licitly 
satisfies boundary conditions (2.8) and H n (±L) = H. Moreover, the y-indepcndcnt quantities H n in ( 2.18| ) have 
clear physical meaning: Being solutions of ( 2.17 ) with = 0, they describe distribution of the local magnetic field 
within the I-layers in the homogeneous Meissner state. (See section III of RefJl.) [Note that in an infinite layered 

N-l +oo 

superconductor H n = 0, and G{n,m) . . . — > (n, m) . . . , where Goo(n,m) are defined by (A17).] By 



comparing Eqs. (2.18) with Eqs. ( 2.10| ) (where /„ = 1), using the property H n = Hn^„ and (All), we establish the 
symmetry relations 



My) = (f>N-n(y), h n (y) = hN-n(y), H n (y) = H N - n (y). 



(2.21) 



Relations ( 2.21| ) n replace the relations (j)n(y) — 4>n+\{y) = <j>(y) and h n (y) = h n+ i(y) = h(y) of the infinite layered 
superconductor .□ 



By the use of Eqs. (|2.18| ), the Gibbs free energy up to first order in r (T) <C 1 can be given the form 



Q(H) = - 



NW 



r(T) 



2H 2 W 



(N-l) 



N-l N-l 



^ £G(n,m) f 

1 ™ i J 



dy 



n—1 m— 1 



d(f) n {y) d(f) m (y) 
dy dy 



N-l " 

+2 5/ 



, • 2 < l ) n{y) 

dy sin 



Af-l 



4^ J2 $ " 



<f>n{L) — 4> n (—L) 

2tt 



(2.22) 



± r. 



fj,~ n + fi- N+n - /i r ' 



,N—r 



-N 



(2.23) 



where H s is the superheating (penetration) field defined in Eq. (A14), and is the total flux carried by a Josephson 
vortex positioned in the rtth I-layer. (The quantities H s and $„ are thoroughly discussed in section IIIB.) The 
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quadratic form of the ele ctrom agnetic energy [the second line on the right-hand side of Eq. (2.22)] can be diagonalized 
with the help of (|Al|), ( [Al^ ): 



[~] — 1 N 1 N 1 ^ 

1 \2 irn(2k + l) . 7rm(2fc + l) f dcj) n (y) d(j) m {y) 
V 2^ A 2Hi L 2^ sm ^ "in - / dy- 



N 



k=0 



N 



N 



dy dy 



(2.24) 



where [u] is the integer part of u. [In deriving ( 2.24 ), we employed the symmetry relations ( 2. 21 ).] Equation ( 2.24 ) 
explicitly shows that physical solutions in a stack with N S-layers are characterized by [-yl Josephson length scales 
A Jfc = A 2fc+ i (fc = 0, 1, . . . , [f ] - 1) from the set ( gig ). 



III. MAJOR RESULTS 



A. Analysis of the equations for the phase differences 



By differentiation with respect to y, integrodifferential equations ( |2.17| ) reduce to a system of TV — 1 ordinary 
nonlinear second-order differential equations 



d 2 My) 
dy 2 



= — [(2 + e 2 ) sm^i(j/)-sin0 2 (j/)] , 



d 2 Mv) _ i 



dy 2 



= — [(2 + e 2 ) sin <f> n (y) - sin 4> n+ x (y) - sin 4> n -i(y)] , 2 < n < N - 2, 



d ^ N d ~ 2 liy) = ^ [(2 + ^ 2 ) sin^^y) - Bin^-ad,)] 



(3.1) 



with bou ndary condi tions (2.9). As shown in section II, physical solutions to (3.1) necessarily obey the symmetry 
relations fl2.li|) and ( pT2l| ). 

A detailed analysis of the analytical structure of the coupled static SG equations ( |3.l|) is the main task of this 
subsection. Such_differential equations were first derived within the framework of a less rigorous, phenctmenological 
approach in Ref.tl . Particular numerical solutions (for H = 0) were obtained in several publications .BtrLl However, 
any analytical investigation of these equations has not been carried out up until now. 



1. Existence and uniqueness of the solution to the initial value problem 



Consider Eqs. (3.1) on the whole axis — oo < y < +oo. Two simple properties of (3.1) are quite obvious: If 4> n (y) 
(1 < n < N — 1) is a solution, the functions 4> n {y) given by 



and 



(f> n (y) = 4> n (y) + 2nk {k is an integer), 



t>n(y) = (frniy + c) (c is an arbitrary constant) 



(3.2) 



(3.3) 



are also solutions. [The latter is a resu lt of the fact that y does not enter explicitly the right-hand side of ( |3.l| ).l Our 
conclusions about the solution to (3.1) will be substantially based on another key property which we formulate as a 
lemma: 



Lemma 1. Consider an arbitrary interval / = \L\^L^ and yo G /. The initial value problem for Eqs. (3.1) with 
arbitrary initial conditions 4> n {yo) = ct n , -^(yo) = Pn has a unique solution in the whole interval /. This solution 
has continuous derivatives with respect to y of arbitrary order and continuously depends on the initial data. (For the 
proof of Lemma 1, see Appendix B.) 

It is worth noting that the existence and uniqueness of a smooth solution to the initial value problem in the whole 
interval / is rather nontrivial for nonlinear differential equations: For such equations, theorems of existence and 



G 



uniqueness are usually valid only locally, in the neighborhood of initial data.Q In our case, global char acte r of the 
solution and its infinite differentiability are ensured by the fact that (f>„ enter the right-hand side of Eqs. (3T) only as 
arguments of the sine. Note that because of the arbitrariness of the interval /, the solution can be uniquely continued 
onto the whole ax is — o o <y < +00 .13 

Differentiating (2.18) with respect to y yields 



N-l 



sin </>„(?/) = e 2 G(n,m)- 



dy 2 



(3.4) 



Multiplying ([TJ) by summing over the layer index n with the use of ( |A1C| ) and performing integration, we 



arrive at the first integral of Eqs. ( p.l 



JV-l 

E 

n=l 



N-l JV-l 



n— 1 m—1 



d<l>n{y) d(f> m (y) 

dy dy 



(3.5) 



C{H) 



2H 2 



N-l 



(N-l)+ ^cos0„(L), 



(3.6) 



where C (H) is the cons tant of integration. Using (3.5), one can eliminate the electromagnetic-energy term from the 
Gibbs free energy ( 2.22 ). .-. 

The existence of the first integral of Eqs. (3.1) was first established in Ref.B. However, the explicit form of the 
matrix elements G(n, m) and the actual value of C (H) was not determined by the authors of Ref.tl. Unfortunately, 
the existence of the first integral for the infinite (N = 00) layered superconductor was not noticed in .any previous 
publications, which partly explains the difficulties with the minimization of the Gibbs free energy£3'E30E3 [Equation 

N-l +00 

( |3.5[ ) holds for N = 00 too, taking account of the substitution G(n,m)... — > Goo(n, m) . . ., where 



m,n— — 00 



Goo(n, 7ti ) are defined by (A17)]. In the case N — 00, the Gibbs free energy fl (H) should be additionally miiiipized 
with respect to the phases tp n , which by the use of the first integral immediately yields the exact solutionEra with 

<My) = <?Wi(y) = <P(y)- 

For our further consideratio n, w e will need an auxiliary lemma: 

Lemma 2. Consider Eqs. (3.1) on the whole axis —00 <y < +00. Conditions 



4>n{ya) = 0, 



dy 



n = 1,2,..., JV-l, 



(3.7) 



where yo ^ ±00, cannot by satisfied simultaneously for all n by any nontrivial solution of (3.1). 

The proof of Lemma 2 is straightforward: Indeed, conditions (3.7) specify the initial value problem for ( |3.l| ). This 
initial value problem has the trivial solution <f>\ = <j>2 = ■ ■ ■ = 4>n-i = 0. By Lemma 1, this solution is unique and 



there are no nontrivial solutions that satisfy (3.7). 

The importance of Lemma 2 lies in the fact that it prohibits the existence of any topological (vortex) solutions 
in any finite interval [L\, L2] in the absence of an external field (H = 0). Indeed, any such solution must necessarily 
satisfy the boundary conditionsuj 4> n (Li) = and (j) n {L2) — 0mod27r for all n, with 0,i(L 2 ) 7^ for at least one 
n = m. For H = 0, we also have ^^(£1) = ^^-(L^) = for all n. However, the former and the latter conditions are 
incompatible by virtue of Lemma 2. 

Another important consequence of Lemma 2 is that for any finite interval [— L, L] the constant of integration in 
(|3.5D satisfies the inequality 



C(H) > N-l, 



(3.8) 



for any H > . Indeed, for the Meissner solution we have 4> n {y) = —<Pn{—y) and </>n(0) = [see (2.11)]. By sett ing 
y = in (3.5), applying Lem ma 2 and using the fact that the qu adra tic form on the right-hand side of ( |3.5| ) is 
positively definite, we get (3.8). For topological solutions, inequality (3.8) is quite obvious. 
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2. Localized solutions. The criterion of existence 



Alongside solutions to fl3.l|) in a finite interval [— L,L], we will consider physical solutions to (3.1) in the infinite 
(—00 < y < +00) and a semiinfinite (0 < y < +00) intervals. To ensure the convergence of the integrals over y 
in the expression for the free energy [see ( 2.22| )1, such solutions must necessarily satisfy the asymptotic boundary 
conditionsc3'Ej 



— ■ > mod 2tt, 

y >±oo 



(3.9) 



d(f> n {y) ^ 

dy y — >±c 



0, 



(3.10) 



for any 1 < n < N — 1. (Note that Lemma 2 does not obtain for y = ±00.) These solutions with a square-integrable 
derivative will be called localized. Their consideration can be reduced to the consideration of the standard initial value 
problem at y = 0, w ith th e ini tial da ta s ubject to a certain existence condition. 

By inserting (3.E) and ( 3.1C ) into ( |3.5| ), we find th e va lue of the constant of integration: C (0) = N — 1. [Compare 
with inequality (3.8) for L < 00, H > 0.] Thus Eq. (3.5) becomes 



N-l 
n=l 



>,<i>n{y) 



e 2 ^ ^ dfin (y) d(j) m (y) 



2 4 ^ ^ -v-'-/ d d 

n=l m=l " y 



(3.11) 



Substituting initial values 0„(O) = a n , ^^-(0) = f3 n into ( ^.ll| ), we obtain the desired condition on a n and f3 n that 
we call the criterion of the existence of localized solutions: 



N-l 



2 N-l N-l 



ri—l m— 1 



(3.12) 



Indeed, Lemma 1 guarantees the existence, uniqueness and differentiability of a solution for arbitrary a n and (3 n in 
the whole interval (—00, +00) [or [0, +00)]. Ow ing to our special choice of the constant of in tegr at ion in ( 3.11 ), the 
solution determined by a n and (3 n obeying (3.12) will necessarily satisfy asymptotic conditions ( |3.9| ), (3.1C). Moreover, 
this solution will automatically satisfy the conditions 



d k (f> n (y) 

dy k 



y — > ±00 



0. 



(3.13) 



f or an y 1 < n < N — 1 and 2 < fc, which can be verified by repeated differentiation of (3.4) and application of (p.Sj), 

Note that loc alized solutions are characterized by the absence of any electromagnetic effects at y — > ±00. Further- 
more, owing to ( 3.11 ), the Josephson energy of localized solutions Ej is exactly equal to their electromagnetic energy 
E em - (From a field-theoretical point of view, this fact is a manifestation of the virial theorem.cj) 



3. Solutions with periodic derivatives. The major Theorem 

Our conclusions about the type and the properties of physical solutions for H > will be based on the major 
Theorem: 

Theorem. Consider Eqs. (fO) on the whole axis —00 <y < +00. The initial value problem 



Myo) = 0, 



dy 



(y ) = 2H>0, n=l,2,...,JV-l, 



(3.14) 



for (3T) has a unique solution on the whole axis — 00 <y < +00. This solution is characterized by the properties 

(t> n (y + P) = (j> n (y) + 2tt, (3.15) 



^( 2/ + P) = 0„( y )>O, 
dy 



(3.16) 
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=2y/H* + H*, foralln=l,2,...,iV-l, 



where if s is defined via ( A14 ), and the period P is given by 

( "2 ' "2 ' ■ ■ ■ ' "2 ) 

P(F) = 



ds 



(0,0,. ..,0) 



N-l 



1 N-i h 2 +h 2 S sin 2 # n 

s n— 1 



(3.17) 
(3.18) 



(3.19) 



ds = eH s 



N-l N-l 



\ n—1 m—1 



The integration in ( |3.19 ) is performed along the integral curve: 



<l>n(y) 



yo<y<yo + ^P, n = 1,2, . . . ,N - 1, 



with (is being the length increment. (See Appendix C for the proof.) 

An immediate physical consequence of the Theorem is the absence of single- vortex solutions to (3.1) in a finite 
interval \—L, L] for any H > 0. (For H = 0, the absence of single-vortex solutions in a finite interval follows from 
Lemma 2.) Indeed, if such solutions existed, they would satisfy at a certain H = Hi > the boundary conditionsE 2 ] 



,(-£) = 0, -p-(-L) = 2Hi>0, n= 1,2,..., N-l, 
ay 



ML) = 2ir, 



dy 



(L) = 2H t > 0; ML) = 0, 



dy 



■(£) = 2Hx > 0, n ^ Z. 



However, these boundar y co nditions are mutually incompatible: The boundary conditions at y = —L specify the 
initial value problem for ( |3.l|) . By the Theorem, this initial value problem admits a unique solution with > on 

the whole axis — oo <y < +oo for all n, whereas the boundary conditions at y = L imply that with n ^ I change 
the sign twice in the interval [-L, L\. 



Analogously, one can prove the absence for H > and, 
well a s such vortex configurations as the "triangular lathee 
(3.15)-(p.l9|) are inherent to the vortex-plane solutions. Ill- 



"l 



oo of any other incahcrcnt vortex solutions, as 
and the "row mode".B In contrast, the properties 



Note that P(H) — > oo, in agreement with Lemma 2. On the other hand, for H 3> H s , the path of integration 

H — >0 

in j3.1Sf) is a straight line: 9 n = J (y - y ), Vo < y < Vo + \P , n = 1, 2, 
For a single thin-layer junction (N 



. , N - 1, which yields P{H) 



In agreement with Rcfs.E 

P{H) = 



equations (3.1) reduce to a single SG equation for 
, the period ( 3.1 9| ) becomes 



j>, with Xjq 



-,K 



(3.20) 



where K (fc 2 ) is the complete elliptic integral of the first kind.EEl Relations (3.15)-( |3~T9 ), anrjlied to cj), determine an 
exact vortex solution in terms of the Jakobi elliptic functions am(«) and dn(u) = £am(«) 



The true vortex-plane solutions appear in a double-junction stack (N = 3). Owing to the symmetry relations (2.21), 
equations (3.1) again reduce to a single SG equation for (pi — (j)2 = 4>, with A,/o = H^T 1 = -7===^. The period P is again 

given by ( 3.20| ) (with redefined H s ), and the exact vortex-plane solution is again expressed via the functions am(it) 
and dn(u). These two exactly-solvable examples (N — 2 and N = 3) not only demonstrate the power and generality 
of the Theorem, but also sh ow , t ba.t the vortex-plane solutions are a natural generalization of the well-known vortex 
solutions in a single junction.ESEaO (See section IV for a more detailed discussion of the cases N = 2, 3.) 
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B. Meissner solutions. The superheating (penetration) fields H s l and H s 

In a finite interval [-L, L], the Meissner solution is characterized by the relations 

cf> n (y) = -<p n {-y), l<n<N-l, 



(3.21) 



resulting from the general symmetry fl2.1l[ ), and obeys the conditions (2.21). Thus, the Meissner boundary value 
problem is completely specified by the boundary conditions 



dy K 



-L) = 2H > 0, <j> n (0) = 0, 1 < n < N - 1. 



(3.22) 



dy 



Up to a certain field H = H s l, the boundary value problem (3.22) has a unique solution with 

interval [— L, L] and — t t < 4> n [— L ) < for all n. The field H s l is determined from the conditions <fi n (±L) = ±7r for 
all n. [As follows from (3J3), (3^), H s l > E s .] Indeed, for H <C H s the Me issne r solution can be obt ai ned in a closed 



analytical form as a linear combination of [y] exponentials with Xjk [see (2.24)]: According to (3.6), (3.8), this case 
allows of linearization. The existence of a unique Meissner solution for higher values of E follows from continuous 
dependence of the solution on the initial data 4> n (-L) = a n and -g 2 *- (-L) = 2H = f3 n . (See Lemma 1.) The 

Meissner solution ceases to exist at H — H s l: For this field, both ^j^{y) and h„(y) are local maxima at y = ±L. 
[See (HI), Q.] 

Physically, the conditions ip n (±L) = ±7r correspond to the vanishing of the surface barrier induced bw Josephson 
currents [j n ,n-i = 0] and the formation of a vortex plane at the side boundaries of the stackB (Thus, at 

H = H s l, the total flux due to the field penetration through the y = ±L interfaces is exactly equal to the flux 
carried by a single vortex plane.) For H > H s l, only topological vortex-plane solutions are possible. However, the 
vortex-plane solutions can exist and be favorable energetically already at fields H < H 8 l- [See the next subsection.] 
For these reasons, H s l should be identified both with the superheating field of the Meissner state and the penetration 
field for the vortex planes. According to the Theorem, at H = H s l, ^^-(0) = 2 \J H* L — for all n, and hence 



H n (0) = JHJ L - El 
implicit equation 



H sL - y/H* L - [G{n, 1) + G(n, N - 1)]. The field H sL itself is determined by the 



2L = P 



(3.23) 



where P is given by ( 3.19| ), with the path of integration being determined by the solution to the boundary value 
problem <f> n (-L) = —n, (f> n (0) = 0. Note that H s l — ► H s for L » H~ x , which will be explained below. 

The consideration of the Meissner effect becomes simpler, when L ^S> A,/ max = Ajo- In this case, the interval [— L, L] 
can be transformed into [0, +oo) by changing the variable y — > y — L and proceeding to the limit L — > oo. In the 
semiinfinite interval [0, +oo), the Meissner solution necessarily satisfies the conditions 4> n {y) — * 0,^%-^- — > 0, 



and we can use the criterion (3.12) that now takes the form 



N 



N-l 



n=l 



■ ^n(O) & 

2 EV 



(3.24) 



where, by (|A14[) and (|A15|) , ( |A16|) 



2a/1 



,N-1\ 



e(A-l)(l + M ^-i) 



(N- 1) N 



E 

fc=0 



X Jk cot 



,7r(2fc + l) 
2N 



(3.25) 
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Physical interpretation of the quantity H s is straightforward. The maximum value of the left-hand side of (3.24) is 
achieved when all <p n {0) — — 7r and is equal to unity, which corresponds to H = H s on the right-hand side. Hence H s 
should be identified with the superheating (penetration) field H soo (H s = H soo ). Note that H s for N < oo is always 
higher than the corre sponding ficl dEm H a = 1 of the infinite (N = oo) layered superconductor. For AT — *• oo, H s — > 1. 

According to (3.24), equations (3.1) can be linearized for H <C H s . In this limit, the explicit Meissner solution in 
the interval [0, +oo) is 



AH L , (2k 
2 , Ajfecot 

fc=0 



N 



2N 



\)tt . (2k 
sin 



1) nn 



N 



exp 



y 

Xjk 



(3.26) 



H n (y) 



2H 

~w 



A ^ cot 



(2k + l)n . (2k + 1) nn 



k=Q 



2N 



■ sm 



N 



■ exp 



V 



(3.27) 



where H n is g iven by (2.20). The distribution of the currents J n and j n ,n-i can be easily obtained from (2.13), (2.14) 
using (|3~26l) , (^271) . 

It i s inst ru ctive to investigate a transition to the infinite layered-superconductor limit N — 1 3> 2 [e _1 ] (e < 1) in 
Eqs. ( 3.26|) , ( ^.27 ). Thus, for I-layers whose layer index n satisfies the condition 



[e- 1 ] <n< N-l - [e- 1 ] , 



(3.28) 



we have H n = H (fi 
we get 



^JV-n) < ^ Under the condition ( |3.28| ), the main contribution to the sums over /c in ( p.26[) , 



Q3.27D comes from < - 1. For such fe and AT > 1, X Jk 



Xj = 1. Therefore, for n in the interval (3.25 



8ff 



+ 00 



■ exp 



HOE 



fc=0 



(-If 

2/c + 1 



= -2H exp (-y) . 



H„(y) « Hexp(-y), 

in complete agreement with the results of Refs.00 . 

Two final remarks would be in order here. The fact that the Meissner solution in a stack with N S-layers is 
characterized by [yl Josephson lenjgths Xjk was not noticed in any previous publications. Moreover, in contrast to 
continuum type-II superconductors ,E3 the local field H n (y) at H = H s cannot be represented as a sum of a p urel y 
Meissner field and a vortex field, because the principle-of superposition does not hold for the nonlinear Eqs. (j3~l[). 
Unfortunately, the latter point was not realized in Ref.E3 concerned with infinite layered superconductors. 



C. Vortex-plane solutions. The lower critical field H c i 



As explained in the Introduction, by a vortex plane we understand a coherent chain of A^ — 1 Josephson vortices 
(one vortex per each I-layer) positioned in a plane parallel to the coordinate axes x, z: see Fig. 2. (In this plane, all 
-J^ and h n are local maxima.) Such solutions in an interval [—L,L] obey the general relations (2.11) with the same 
constant 0mod27r ^ for all n and the symmetry conditions ( p. 21 ): hence the same set of [-y] Josephson lengths 



Xjk as in the Meissner case. The Theorem envisages the existence of vortex-plane solutions for any H > 0. Using 
the Theorem and Lemma 1, we can easily establish all basic properties of the vortex-plane solutions. They can be 
summarized as follows. 

The boundary value problem for the vortex-plane solutions is specified by the boundary conditions 



~dy- { ~ L) 



2H > 0, 



<t> n (0) - nN v 



1< n < N - 1, 



1,2, 



(3.29) 



where A^ is the number of vortex planes. As in the case of the Meissner boundary value problem ( [3.22 ), the boundary 
value problem (3.2E) for a fixed N v and any L > has a unique solution, with > in the whole interval [-L, L] 
and — 7r < <j) n (-L) < for all n, in a certain field range 
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Hi 



H? < H < H 



(3.30) 



where H Q = H s l. The lower bound of the existence of the N v -vortex-plane solution is determined by the boundary 
value problem 



= 0, K (0) = nN v , l<n<N-l, 
whereas the upper bound is determined by the boundary value problem 

= -7r, (j>n(P) = *N v , l<n<N-l. 
The field -ff/v„ is determined by the implicit equation 



2l = {n v + i)p(Jh% v -h^ 



(3.31) 
(3.32) 

(3.33) 



where P is given by ( 3.19| ), with the path of integration being determined by the solution to the boundary value 
problem (3.32). Note that for N v = (the Meissner state) Eq. ( ft.33j ) coincides with Eq. (3.23), as it should. 
Physically, P deter mines the period of the distribution of H n (y) and j n . n -i{y)- 

As follows from ( 3.30 ), the range of the existence of the state with N v vortex planes overlaps with that of the 
state with N v — 1 vortex planes. Thus, the state with N v = 1 overlaps with the Meissner state. For L ^> if" 1 , the 
overlapping can occur for several neighboring states. The overlapping is negligibly small only for L <C H^ 1 (with 
arbitrary N v ) and for N v 3> 1 (with arbitrary L) . (For ordinary Josephson vortices in a single junction, the overlapping 
is well known r 4 lr 5 l ) However, the actual number of vortex planes, N v , for given H should be determined from the 
requirement that the free energy be an absolute minimum. Note that for L — > oo, yj Hg — i?f = \J H^ L — iJf — > 
and the state with N v — 1 appears already in zero external field. This case is discussed in more detail below. 

For L 3> Hj 1 , we can proceed to the limit L —> oo and consider the state with N v = 1 in the interval (— oo, +oo). 
With the asymptotic boundary conditions 



0, 



M -? 27T, 

y — >-+oo 



d(f> n (y) 

dy 



y — ^ioo 



(3.34) 



for all n, this state satisfies the criterion (3.12) and can be obtained as the solution to an initial value problem at 
y = 0. Indeed, it can be constructed from the Meissner solution in the interval [0, +oo) at H = H s : Using the property 
(3.2), by means of the transformation 0„ — > cf> n + 2n we obtain in the interval [0, +00) a solution that satisfies the 
initial conditions 



dy 



(0) = 2H S , l<n<N-l, 



(3.35) 



and the conditions (3.34) for y — > +00. By Lemma 1, this solutions can be uniquely continued into the whole interval 
(—00, +00). The solution thus obtained satisfies the conditions (3.34) for y — > ±00 and hence represents the desired 
singe-vortex-plane solution. By the construction, h n (0) = H s [1 — G(n, 1) — G(n,N — 1)]. The total flux carried by 
the vortex plane is 



N-l 



J2 = t (N - 1) 



n=l 



2V1 + T- £ 1 



N-l 



e(N-l) l + fx 



N 



(3.36) 



Shy ( 2.23| ). Note that in contrast to Josephson junctions with thick electrodes!^ and infinite layered 
□ the flux carried by a Josephson vortex in a finite thin-layer S /I structure is not quantized and is 



where $„ is gives 
superconductors,! 

always smaller than the flux quantum $ = ft- For N — 1 3> 2 [e _1 ] , $ — > ?r (N — 1), as it should. 

To determine the thermodynamic lower critical field H c \ at which the vortex-plane solutions for L 3> H~ l become 
energetically favorable, we must calculate the difference between the Gibbs free energy in the presence of a single 
vortex plane, Q V (H), and t he Gibbs free energy of the homogeneous Meissner state, SIm(H) [the sum of the phase- 
independent terms in (2.22)]: 



n v (H) - n M (H) 
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r(T) 



N-1N-1 +?° 

G(n,m) / cfy 



(y) #m(y) 

dy 



4$iJ 



(3.37) 



where the total flux $ is given by ( 3.36| ). The first term on the right-hand side of (3.37) should be interpreted as the 
self-energy of the vortex plane: E v = 2E em — 2Ej. From (3.37), we get: 



Hci — 



E v 



4r(T)$ 



(3.38) 



For thermodynamically stable solutions, we must necess aril y have H c \ < H s . It is straightforward to verify that this 
condition is met by the vortex-plane solutions. Using (3^), ( |3.34 ), the initial values (3 n = 2H a and integrating by 
parts, we convert E v into the form 



E v = 2r(T) 



N-l " 

/ dy(j> n (y)Bm<l> n (y)-2H s $ 

™ =1 +oo 



(3.39) 



The first term on the right-h and side of ( 3.39 ) is positive, because in the region < y < +oo all </>„ satisfy the relation 
7T < 4> n < 2tt. By the use of (3^) and ( |A14| ), we obtain the following strict inequalities: 



N-l " 

2H S <S> < Y J d V<t>n{y) sin 4> n (y) < 4iJ s $, 



Hence, 



< E v < 4r(T)H s $. 



< H cl < H s , 



as anticipated. Note that in all special cases admitting exact analytical solutions (N = oo,Ela and N = 2,3, see section 
IV), H cl = \H S . 

D. Single-vortex solutions for H = 0, L — oo, and other localized incoherent vortex solutions 



As shown in section IIIA, the only topological (vortex) solutions admitted by Eqs. (pTj) for H > are the vortex- 



plane solutions. For H = 0, equations ( 3.1 



Therefore, incoherent vortex solutions to (3.1 



do not possess any topological solutions in a finite interval [— L, L]. 
L) can exist only for H = 0, in the infinite interval (— oo, +oo), and must 
necessarily meet the requirements imposed on localized solutions. Below we consider in detail the most important 
type of such solutions, namely single- vortex solutions. 



A single Josephson vortex positioned in the Ith I-layer at y = obeys the symmetry relations 

<j>i(y) = 2tt- <t>i(-y); (f> n (y) = -<t> n {-y), n^l, 



see (2.11)] and the asymptotic boundary conditions 



k{y) -» 0, Mv) -? 2tt; 

y — > — oo y— »+oo 



y— »±oo 



(3.40) 
(3.41) 



d<t>n{y) 

dy y^±c 



0, for all 1< n < N - 1. 



(3.42) 



Moreover, dh ^ v ^ > in the region — oo < y < 0, and dh 2y V ^ < in the region < y < +oo. Hence, <p n must satisfy 



the relations 



dy 



< <f> n {y) < 7T, y e (-oo, 0) ; -ir < <j> n {y) < 0, ye (0, +oo) , for n ^ I, 



(3.43) 
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and the initial conditions 



ai = <pi(0) = it: 



</>„(0) = 0, n^l, 



(3.44) 



_ #K0) #»(0) 

Pi = —j > 0; f) n = — 

ay dy 



< 0, I. 



(3.45) 



A necessary condition of the existence of such a solution is provided by the general criterion (3.12) and has the form 

2 N-l N-l 



(3.46) 



The flux through the nth I-layer due to the vortex in the Ith I-layer can be found using and ( |3~4l|) 



+00 

®ni = / dyhniy) 



TIC 



2a/1 



|n-H| _ ^ I/ 1 M j + A* (M "J 



The total flux carried by the vortex in the Ztli layer is $; = Yl n =i ^nl, where $; is given by ( 2.23 ) with n = I. For 
N < oo, the total flux $/ is not quantized and is less than the flux quantum $ = (See the previous subsection.) 

In contrast to the Meissner and the vortex-plane solutions, single-vortex solutions, in general, do not obey the 
symmetry (2.21), inherent to the original integrodifferential equations (2.7), (2.17) minimizing the Gibbs free-energy. 
Therefore, they are characterized by the full set ( A15 ) of N — 1 length scales A^. (The only exclusion is a stack with 
an odd number of junctions N — 1 and I = y.) 

Although the energy of single- vortex solutions is lower than the self-energy of the vortex-plane solutions at H = 0, 
the former solutions cannot be used, e.g., for estimates of H c \ because of their absolute instability at H > 0. 
Unfortunately, this crucial issue was not understood in any previous publications. Thus, the main distinctive feature 



(3.47) 



of single-vortex solutions that precludes their existence at H > 0, i.e 



Refs 



< for n ^ I, was overlooked in 



concerned with infinite layered superconductors. 



points of view in Refs. 



.) The negative sign of 



rf0n(O) 



(The approach of Refs. was criticized from different 

.111 . However, the effect of this 



, for n ^ I was noticed in Ref.t 
property on the stability of single-vortex solutions at H > was nat_cealized therein. 

The feature ( 3.45| ) is clearly reproduc ed b y the numerical studigasttffl of finite stacks. Unfortunately, in the absence 
of any analytical investigations of Eqs. (3.1), the authors of R efejMp Lj failed to establish the actual domain of validity 
of their numerical results: H = 0, L = oo. The condition (|3.46| ) was not found either. In this regard, it would 
be appropriate to point out an intrinsic limitation of the numerical calculations that does not allow one to consider 
them as a proof of the existence of single-vortex solutions even in the case H = 0. By necessity, these calculations 
are performed in a finite interval [Li,^]- However- q« shown in section IIIA, equations (3.1) do not admit any 
topological solutions (in a strict mathematical senseoo) in any finite interval at H = 0. It should be also noted that 



the single- vortex solution cannot be represented by a linear combination of the terms 4 arctan exp 
the principle of superposition does not apply to the nonlinear Eqs. 



(3.1) 



because 



The consideration of other localized incoherent vortex solutions to (3.1) (i.e., solutions with 2 < k < N — 1 vortices 
in the plane y — 0, and vortex-antivortex configurations) can be done along the same lines. At H > 0, all these 
solutions are unstable too. 



IV. THE TWO EXACTLY-SOLVABLE CASES 



Below, we present an exact solution for the cases N — 2,3. Valid for any H > and L > 0, this solution provides 
a clear illustration of the general results of sections IIIA-IIIC. The single-junction case (N — 2) illuminates some 
typical specific features of the thin-S-layer limit and establishes a profound physical and mathematical relationship 
between the ordinary Josephson vortices and the vortex-plane solutions. In the double-junction case (N = 3), our 
exact solution describes the true vortex planes, which is of considerable experimental interest: For a double-junction 



stack, we already have direct o 
with the vortex-plane solutions. 



aervations of Joscphson-vortex configurations that can be unambiguously identified 
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A. A single thin-layer junction (N = 2) 



In this simplest case, a single phase difference <j)\(y) = <fi(y) satisfies the usual static SG equation 

dV(v) _ i 



dy 2 



A 



2 sin0(y), 

,70 



with the Josephson lengtl 



>,70 



(4.1) 



(4.2) 



Note that A,/o, given by (4J2), for e <C 1 is much smaller than the J oseph son le ngth of a single junction with thick 
electrodes, which in our dimensionless units is Aj = -^/^.Ell From ( 2.20| ) and (3.25), we get the local field in the 
homogeneous Meissner state 



2H 
2+1 2 



and the superheating (penetration) field 



H » - A J0 = 



V2- 



(4.3) 



(4.4) 



respectively. For e <C 1, the superheating (penetration) field H s , given by (4.4), is much higher than the corresponding 
fieldBO of a single junction with thick electrodes H s — A,/. 



Using the Theorem, we immediately find an exact solution to the boundary value problems ( 3.22 ) and ( 3.29 ) 

y 



(j)(y) = 7T (N v - 1) + 2am 



kX 



K(k z 



.70 



(4.5) 



dn 



A' A 



.70 



kXjoH 



N„ 



2m, m = 0, 1, . 



(4.6) 



0(y) = + 2am 



( V 



10 



(4.7) 



dn 



fcA 



70 



k\ Ja H, 



N v = 2m+ 1, 



0,1,. 



(4.8) 



where_am(7t) and dn(u) = £am(«) are the Jakobi elliptic functions, and K (A: 2 ) is the elliptic integral of the first 
kindB The range of the existence of the solution with N v = 0, 1, . . . Josephson vortices is given by ( |3.30[ ), where Hjy^ 
is determined by the implicit equation 



-K 



A jo H ff v A jo V ^iv„ ^ jo 
Note that although the Meiss ner a nd t he vortex solutions to (PO) were studied a long time ago 



(4.9) 



closed-form 



analytical solution of the type (4.5)-(4.9), valid for any < L in the whole field rapSg < H <C H c2 , has not been 
found up until now. (Apparently for this reason, there exists an erroneous beliefM] that Josephson vortices "do 
not fo rm" in single junctions with small L.) Using asymptotics of am(u), dn(u) and K (fc 2 )JH3 we can obtain from 
(4.5)-(4.9) all physically interesting limiting cases. 
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1. The Meissner solution for [0 < y < +00) 



For the fields < H < H s 



+e , the Meissner solution in the semiinfinite interval [0, +00) can be obtained from 
(4.5), (4.6) with N v = by changing the variable y —* y — L and proceeding to the limit L — > 00. Using the formulas 
of section II, up to first order in r(T) <C 1, we have 



4>(y) — — 4arctan - 



H exp 



H(y) = H 1 (y) = h(y)+H u 



(4.10) 
(4.11) 



h{y) = hx{y) 



2X JQ H H s + y/H* - H 2 exp 



H s + yJE* - H' 2 + H 2 exp 



A, 70 



(4.12) 



j(y) = h,o(y) = H s + ^H 2 - H 2 



[H s + V-fff - H 2 


2 

- H 2 exp 




exp 


y_ 

A, 70 




H s + ^H 2 - H 2 


2 

+ H 2 exp 


" ?y_~ 

A, 70 


1 2 



(4.13) 



J(y) = J (y) = J^y) = —[H-H x - h(y)} 



f(y) = f (y) = h(y) = l- 



r(T) 



^h{y) + ^[h{y) + H x -HY 



(4.14) 
(4.15) 



2. The single-vortex solution for (—00 < y < +00) 

The single-vortex solution in the infinite interval (— 00, +00) can be obtained from fl4.7| ), (| 
proceeding to the limit L — + 00. It has the form 



<p(y) = 4arctanexp 



>J0 



with N v — 1 by 
(4.16) 



h(y) = ^JO cosh 1 



V/o 



(4.17) 



j(y) = -2 cosh" 



V70 



sinh 



.70 



(4.18) 



The quantities H(y) , J(y) and f (y) are given by ( 4.11 ), ( 4.14 ) and ( 4.15 ), respectively, with h{y) taken from ( 4.17 ) 
By inserting (4.16) into (3.37) with H = 0, we obtain the vortex self-energy: 



E v = 8r(T) 



(4.19) 
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The vortex flux, according to ( |3.36| ) , is 



$ = 7T- 



2 + e 2 



and the lower critical field, by ( 3.38 ), is 



Hd — —H s — 



2 + ? 



(4.20) 



Thus, for e <C 1, the vortex flux $ <C $o = an d the lower critical field (4.20) is much larger than the corresponding 
field of a single junction with thick electrodes H c \ = 



in agreement with Ref.tj. 



3. The vortex solution for small Josephson screening 

When the screening by Josephson curren ts is negl igibly small, i.e., (i) for L -C Ajo and arbitrary H, or (ii) for 
H a <C H <C -ff C 2 and arbitrary L, equations ( |4.5| )-( |i~9| ) become 

(f>(y) = 7riV„ + 2i?y 



(-1) 



[sin (2i?y) - 2Hy cos (i/W 7 )] 



(4.21) 



where N v = [^] . Equations of the type ( |4.2lD were first obtained for the infinite layered superconductor by means 
of a perturbation theoryBcl As can be seen from (4.21), the overlapping of states with different N v now vanishes, and 
the period of the vortex structure is given by P — -g>, in full agreement with the general results of sections IIIA and 
IIIC. 



B. A double-junction stack (JV = 3) 



Owing to the symmetry (2.21), 4>%(y) — 4>2{y) = 4>{y), an d a double-junction stack is described by the single SG 
equation ( p~l| ) with the sing i 



c Josephson length 



VT + ? 



(4.22) 



An exact solution to the boundary value problems ( |3.22[ ) and ( [3.29 ) is again given by Eqs. (4.5)-(4.9), where now N v 
should be interpreted as the number of vortex planes. Thus, practically all the formulas of the previous subsection 
(with corresponding redefinition of the physical quantities) hold for the double-junction stack too. Note that 

Mv) = My) = <t>{v), fli(v) = H 2 ( y ) = H( y ), h 1 { v )^h 2 {y) = h( y ), 
hAv) = h.i(y) = j(y), My) = My) = J {y)-, My) = My) = f{y), 



r(T) H 

My) = o, My) = i - -rr-Hy), h 1 = h 2 = 



vo 



1 + e 



2 ' 



According to ( 3.25| ), the superheating (penetration ) field is 



H„ = A 



70 



which is smaller than the corresponding single-junction value (4.4), in agreement with Ref.0. 
The vortex-plane self-energy for the interval (— oo, +oo) is 
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E v = 16r(T)A J0 = 16r(T) 



VTT- 



(4.23) 



and the flux is 



$ = 2tt 



1 + e 2 ' 



which immediately leads to the lower critical field: 



Hd — —H s 



2 vT+e* 



As can be seen by comparing ( 4.23 ) with the single-junction expression ( P~19l) , the energy per vortex in the double- 
junction stack is higher. The vortex-plane solutions are presented schematically in Fig. 2. 



V. DISCUSSION 



Within the framework of standard methods of the theory of differential equations, we have obtained a complete 
mathematical description of the Meissner effect and the vortex structure in periodic thi n- lay er (N — Injunction stacks 
(2 < N < oo). The results of our analytical analysis of the coupled static SG equations (3.1), summarized in Lemmas 
1,2 and the Theorem, should provide a basis for any further analytical or numerical study of these equations, not 
necessarily restricted to the field of superconductivity. 

By proving the absence of single- vortex solutions to (3.1) for H > and specifying .the actual domain of their 
existence (H = 0, L = oo), we have clarified a wide-spread misunderstandingOBlIEEj The previous estimates 
of H c i turn out to be irrelevant because of absolute thermodynamic instability of single-vortex configurations. In 
full agreement with our study of infinite (N = oo) layered superconductors,™! we have shown that the only true 
topological (vortex) configurations that survive at H > are the vortex planes. Being a natural generalization of 
ordinary Josephson vortices in a single junction (iV = 2), t he vortex-plane solutions for 3 < N inherit such properties 
of the former as period icity along the layers [Eq. 



(3. If 



and the overlapping of states with different topological 
numbers N v [Eq. ( 3.30| )]. A unified mathematical approach of this paper, valid for 2 < N < oo, allowed us to derive 
an exact expression [Eq. ( 3.25| )] for the superheating (penetration) field H s as an explicit function of N. Within 
the framework of the same unified approach, we have obtained a new exact solution to the single SG equation, Eqs. 
( 4.5 )-( 4.9 ) , valid for any H > and L > 0. This solution refutes the assertionsHQ that Josephson vortices "do not 
form" in single junctions with small L. 



Our closed- form analytical expression for the Meissner field, Eq. (|3.27| ), clearly illustrates an important feature of 
the Meissner effect in (N — Injunction stacks, not noticed in previous theoretical publications, namely the existence 
of [y] different Josephson -Lengths Xjk (k = 0, 1, ... , [y] — !)■ This result may prove to be usefuiin view of the 
current experimental effortsciriEl to verify the interlayer tunneling model of high-T c superconductivitytia by measuring 
the c-axis penetration jdepth. [The penetration of the parallel magnetic field with a distribution of length scales has 
been recently observedEa in the organic layered superconductor k-(BEDT-TTF)2Cu(NCS)2-] 

Finally, the vortex-plane solutions of our paper provide rSunatural explanation for the observed coherent Josephson- 
vortex configurations in artificial double-junction stacksll3 and weakly-coupled multilayersliia in the presence of a 
static parallel external field H > 0. The observation, of isolated interlayer vortices in the high-T c superconductor 
Tl 2 Ba2Cu0 6+<5 y and in K-(BEDT-TTF) 2 Cu(NCS) 2 Eil can be explained by the fact that in these experiments the 
external field was set equal to zero (H — 0). In this situation, isolated Josephson vortices can be pinned by structural 
defects, because their self-energy is lower and the spatial extension along the c-axis is smaller than those of the vortex 
planes at H = 0. We hope that our results will stimulate further experimental investigations. 



APPENDIX A: THE SOLUTION OF THE FINITE DIFFERENCE EQUATION FOR Hn(Y) 



Equations ( [2.17 ) can be regarded as a nonhomogeneous finite difference equation for H n (y) witi respect to the 

According to general theory of such equations,E3 its solution can 



layer index n, subject to boundary conditions (2 
be represented in the form 



H n (y) = h n (y) + H„, 



(Al) 



where 
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JV-1 



K{y) = y XI G ( n > '' 



711 — 1 



dy 



is the particular solution of (2.17) satisfying the boundary conditions 

ho(y) = h N (y) = 0, 

and 



(A2) 



(A3) 



(A4) 



e 2 / e 2 

M = i + -- e yi + -, (o< M <i), 



(A5) 



is the solution of the homogeneous form of ( 2.17 ) (with the zero right-hand side) meeting the boundary conditions 

H a = H N = H. (A6) 



The quantities G(n,m) in (A2) are elements of a (A — 1) X (N — 1) matrix G(l < n,m < N — 1). They obey the 
nonhomogeneous finite difference equation 



(2 + e 2 ) G(n, to) - G{n + 1, to) - G(n - 1, to) = 5 n , r 
(S n ^ rn is the Kronecker index) with the boundary conditions 

G(0,to) = G(A,m) = 0. 

The explicit form of G(n, to) is 



G(n, to) 



2eJl 



The following properties of G(n,m) can be easily verified using ( |A7| ) and (|A9|): 

G(n, to) = G(m, n), 



(A7) 
(A8) 

(A9) 
(A10) 



G(n, N -m) = G(N - n, to), 



G(n, m) > for any 1 < n, m< N — 1, 



(All) 
(A12) 



iV-l 

X G(n, to) = - [1 - G(n, 1) - G(n, A - 1)] 



1 < n < N - 1, 



(A13) 



iV-l N-l 



m— 1 



A- 1 - 



1 + T 



,AT-1 



A' 



A- 1 



(A14) 



According to (A7), the inverse of the matrix G(n, m) is a Jacobian matrixo G 1 (n, to) with elements G 1 (rt, n) = 
2 + e 2 (n = 1, 2, . . . , A — 1), G^n + l,n) = G- 1 (n,n + 1) = -1 (n = 1, 2, . . . , N - 2), and G^^m) = for 
|n — to | > 1. Hence all the eigenvalues ej of G(n, m) are positive and given by 
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A 2 

ei = 4> X 3 = —< i = l,2,...,N-l, (A15) 

■h + e 2 - 2 cos ZL 



with the corresponding eigenvectors 

2 . Trij 



"'^ViV^T' i,i = l,2,...,iV-l. (A16) 

The quantities Aj determine the characteristic length scales of the stack with N S-layers. 

Note that for an infinite layered superconductor, the analog of the matrix G(n, m) is Goo(rt, m), determined by the 
matrix elements 

Jn-m| 

Goo(n,m) = , — oo<n,m<+oo, (A17) 



that satisfy the summation rule 

JV-l 



^ Goo(n, to) = i. 



APPENDIX B: PROOF OF LEMMA 1 

By introducing new functions 

ipiiv) = <l>i(.v),ip2(y) = <fa(y), ■ ■ .,ipN-i{y) = </>n-i(v), 
, t \ d4> x {y) d<p 2 {y) , , s d<f>N-i(y) 

tpN(y) = — -, — ,ipN+i{y) = —j — , • ■ ■ , ip2N-2(y) = — -, , (Bi) 

ay dy dy 



we convert (3.1) into an equivalent normal system of 2N — 2 first-order equations 

— ^r^- = Fi (ipi,ip2, i>2N~2) , l<i<2A^-2, (B2) 
dy 

Fi (01,02, ■ ■ - ,1p2N-2) = V'i+JV-lj 1 <i < N -1, 

Fn (ipi,ip2, ■ ■ ■ ,ip2N-2) = -a [( 2 + e 2 ) sin -0i - sin 02 ] , 

F; (i/>i,ip 2 ,'--,ip2N-2) = [( 2 + e 2 ) sin^i -sin^-i -sin^+i] , N + 1 < i < N - 3, 

F 2 jv-2 ("01, "02, ■ ■ .,ip2N-2) = -2 [( 2 + e2 ) sin0Ar„i - sinV'jv-2] , 
subject to initial conditions 

0i(2/o)=o;i, l<*<iV-l, 

iMl/o) = A-JV+i, Ar<z<2Ar-2. (B3) 

To prove the statement of Lemma 1, it is sufficient to observe that all Fi (0i, 02) ••• , 4>2N-2) arc continuous 
functions of their arguments for y S (—00, +00) and 0^ £ (—00, +00) (1 < k < 2N — 2). Moreover, their partial 
derivatives with respect to 0^ satisfy the relation 
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for y E (— oo, +00) and ipk € (—00, +00) (1 < i, k < 2N — 2). Thus, the Lipschitz conditions with_respect to ipk 
are met for y S (—00, +00) and tfj^ £ (—00, +00) (1 < k < 2N — 2), which immediately guarantees^ the existence 
and uniqueness of a solution to (B2), satisfying arbitrary initial conditions (B3), in an arbitrary interval / = [L±, L2] 
such that yo S I. Continuous dependence of the solution on initial data is a result of continuous dependence 
of Fj, (ip±, ip2i ■ ■ ■ i ^2N-2) on their arguments and of the condition (B4). Infinite differentiability of the solution 
automatically follows from infinite differentiability of Fi {ipi,ip2, ■ ■ ■ ,4 } 2N~2) with respect to their arguments. 



APPENDIX C: PROOF OF THE THEOREM 



Owing to the property (3.3), without loss of generality, we can consider the initial value problem (3.14) for y — 
ya = 0. To p rove the Theorem, we have to show that there exists a solution to ( |3.1| ) in (—00, +00) with the properties 
( 3.15 )-( [3.19| ) that at y = yo = s atisfies ( 3.14 ). By Lemmal, this solution will represent the sought unique solution 
to the initial value problem (3.14). 

Consider an arbitrary finite interval I = [— ^ , ^1 . We start with the Meissner boundary value problem in the 
interval /: 



dy 



±| ] =2/7 -.(. 



MO) = 0, 



1< n< N - 1. 



(CI) 



Up to a certain H = H s , the problem (|C1|) has a unique solution with the properties 

My) = -M-v)> y^ 1 ' i<n<jv-i, 



(C2) 



-P(y)>0, ye I, l<n<N-l, 
dy 



(C3) 



7T < 



^)<0, O<0„[-)<^, l<n<JV-l, 



(C4) 



where H s is determined by the conditions 



) = ±7r > 



1< n < N- 1, 



(C5) 



and satisfies the inequality H s > H s , with H s defined via (A14). (See section IIIB.) 

Now we will prove that all ^p- (0), where 4> n are the solution to (plf), for H — > H s tend to the same limiting value 



dy 



(0) = 2<JHj - HI \<n<N-l. 



(C6) 



Indeed, relations (|C6]) satisfy (3.5) at y — with (3.6), where the substitution H — > _ff s , L — > has been made. 



However, it is still necessary to show that the limiting value (C6) is uniquely determined by the solution to (|Cl|). To 
this end, we consider the relation 



s n=l m=l n—1 



N-l N-l 



d(3 n 



N-l 



P\ #„ (f ) 



> 0, 



(C7) 



where /3„ = ^ L (0). [Relation flC6| ) is obtained from (3.5), ( ft .6] ) by the substitution — > H , L — > and differentiation 
with respect to H.j Notice that 



/3 n <2H, l<n<N-l, 



(C8) 
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in the whole range < H < H s , because all have minima at y = 0. [See (3.4).] With increasing H from zero to 
H s , all <p n (-j) increase monotonously from zero to 7r. Thus, the difference on the left-hand side of ( |C6| ) first increases 
(for < H <C H s ) and then decreases (for H s < H < H s ). It means that the dependence f3 n (H) changes from a 
linear one for < H <C H s [i.e., /3 n (H) = b n H, < b n < 2, 1 < n < N — 1] to a nonlinear one for H s < H < H s , with 



^>2, 
dH 



1< n< N - 1. 



At H = H B the right-hand side of (pq) is equal to zero, and we have 



where 6„ satisfy the conditions 



N-l N-l 
77 = 1 m—1 



G(n, m)b n b ri 



4(JV- 1) 
^1 : 



(C9) 



(CIO) 



(en) 



< 6 n 



< 6„ < 6 n 



=min {6„} < 2, 



=max {b n } > 2. 



The integration of ( CIO ) for n = m under the condition [3 n (H s ) — yields 



With the help of ( |C8| ) and (|C9| ), taken at H — H s , we establish the following inequalities for 6 m j n and 6 D 

2 



i? 2 

1 -«!<"»»' 



< 2, 



2 < 6„ 



< 



H 2 



(C12) 



(C13) 



(C14) 



H 

Inequalities (C14) must hold for any H s > H s . Proceeding to the limit — f- — * 0, we get & min = 6 max = 2, hence the 



Hi 



relations (C6). 

The Meissner solution at H = H s can be uniquely continued from the interval / into the whole interval (— oo, +oo).EI 
The solution (j> n (1 < n < N — 1) thus obtained possesses all the required properties (3.14)-(3.19), where yo = 0, and 



H = y Hi — iff. To show this, it is sufficient to prove the property ( p.lE| ). [The properties ( 3.14 ), ( 3.17 ) and ( 3. IS ) 
are o bvious. The property (3.16) results from the differentiation of ( 3.15| ). The property ( 3.19| ) is obtained from (|3.5|) . 
( 3.6) by integration.] 

Consider a set of functions 

My)=^n{y + P)-2ir, l<n<N-l, (C15) 
where 4> n {y) (1 < n < N — 1) are the continuation of the Meissner solution at H = H s into the interval (— oo, +oo). 



Owing to the properties (|3.2|), (3.3), the functions (C15) also represent a solution to (3.1) in the interval (— oo, +oo 



At y = — -j, the solution (f> n (1 < n < N — 1) satisfies the same initial conditions as the solution <f> n (1 < n < N — 1). 
By Lemma 1, the latter means that <fi n = <\> n (1 < n < N — 1) in the interval (— oo, +oo), i.e., 

cp n {y + P)-2n = 4> n (y), 1 < n < N - I, (C16) 

which accomplishes the proof of the Theorem. No te that in this proof we have not employed the symmetry relations 
( [2.21 ) resulting from the boundary conditions 
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FIGURE CAPTIONS 



1. Fig. 1. The geometry of the problem (schematically). Here N — 4, a <C p (a is the S-layer thickness), 2L = W, 
H>0. 

2. Fig. 2. The vortex-plane solutions in a double-junction stack (schematically): a) the state with a single vortex 
plane (N v = 1); b) the state with two vortex planes (N v = 2). The black dots conventionally denote the 
positions of Josephson vortices, and the arrows show the distribution of intralayer and Josephson currents. 
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